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A  linear  analysis  of  the  shallow  water  equations  in  spherical  coor¬ 
dinates  for  the  Turkel-Zwas  (T-Z)  explicit  large  time-step  scheme 
is  presented.  This  paper  complements  the  results  of  Schoenstadt, 
Neta  and  Navon,  and  others  in  1-D,  and  of  Neta  and  DeVito  in  2-D, 
but  applied  to  the  spherical  coordinate  case  of  the  T-Z  scheme. 
This  coordinate  system  is  more  realistic  in  meteorology  and  more 
complicated  to  analyze,  since  the  coefficients  are  no  longer  constant. 
The  analysis  suggests  that  the  T-Z  scheme  must  be  staggered  in 
a  certain  way  in  order  to  get  eigenvalues  and  eigenfunctions  ap¬ 
proaching  those  of  the  continuous  case.  The  importance  of  such 
an  analysis  is  the  fact  that  it  is  also  valid  for  nonconstant  coefficients 
and  thereby  applicable  to  any  numerical  scheme.  Numerical  experi¬ 
ments  comparing  the  original  (unstaggered)  and  staggered  versions 
of  the  T-Z  scheme  are  presented.  These  experiments  corroborate 
the  analysis  by  showing  the  improvements  in  accuracy  gained  by 
staggering  the  Turkel-Zwas  scheme.  ©  1997  Academic  Press 


1.  INTRODUCTION 

The  transfer  function  was  first  introduced  in  electrical 
engineering  by  Stremler  [7].  It  is  different  from  Fourier 
analysis  since  it  gives  information  about  amplitude  distor¬ 
tion  and  not  just  phase.  Schoenstadt  [2]  has  applied  the  idea 
to  comparison  of  various  schemes  for  the  one-dimensional 
shallow  water  system.  Neta  and  Navon  [3]  extended  the 
results  to  include  the  Turkel-Zwas  [1]  explicit  large  time- 
step  scheme  on  a  limited-area  domain  for  the  shallow  water 
equations.  Steppeler  [8,  9]  analyzed  some  hybrid  methods. 
Neta  and  DeVito  [4]  extended  Schoenstadt’s  results  to  the 
two-dimensional  case.  Neta  [10]  has  extended  the  transfer 
function  analysis  to  the  two-dimensional  Turkel-Zwas  ex¬ 
plicit  large  time-step  finite  difference  scheme  in  a  Cartesian 
coordinate  system.  Song  and  Tang  [13]  have  used  the  La¬ 
place  transform  instead  of  the  Fourier  transform  to  analyze 
the  unstaggered  as  well  as  the  staggered  Turkel-Zwas 
scheme. 

1  Part  of  this  research  was  conducted  while  the  author  was  visiting  the 
Department  of  Mathematics  at  the  Technion,  Haifa,  Israel. 

2  To  whom  correspondence  should  be  addressed. 

3  Current  address:  Naval  Research  Laboratory,  Monterey,  California 
93940. 


In  this  paper  we  extend  the  linear  transfer  function  anal¬ 
ysis  to  the  two-dimensional  shallow  water  equations  in 
spherical  coordinates  for  the  Turkel-Zwas  discretization. 
Actually,  we  show  how  to  obtain  the  modal  expansion  for 
the  shallow  water  equations  in  spherical  coordinates  and 
for  the  Turkel-Zwas  discretization  of  these  equations.  At 
this  point,  we  should  comment  on  the  choice  of  the  method. 
Computationally  efficient  and  accurate  schemes  for  the 
numerical  solution  of  the  shallow  water  equations  are  of 
crucial  importance  in  atmospheric  and  oceanographic 
models.  Two  different  approaches  have  been  taken,  both 
of  them  dealing  with  the  different  time  scales  of  advective 
(Rossby)  waves  and  gravity  inertia  waves  separately.  The 
first  of  these  was  the  split-explicit  and  semi-implicit 
schemes.  The  Turkel-Zwas  scheme  takes  a  different  view 
by  proposing  a  space  (rather  than  time)  splitting  approach. 
This  is  based  on  the  fact  that  the  fast  gravity  inertia  waves 
which  restrict  the  time  step  contain  only  a  small  fraction 
of  the  total  available  energy  and  therefore  can  be  calcu¬ 
lated  with  a  lower  accuracy  on  a  coarse  grid.  The  Rossby 
waves  which  contain  most  of  the  energy  are  calculated  on 
the  finer  mesh.  When  the  ratio  of  the  coarse  and  fine  grids 
is  an  integer  p  >  1,  one  can  use  time  steps  nearly  p  times 
larger  than  those  allowed  by  the  usual  explicit  scheme. 
The  CFL  condition  for  the  Turkel-Zwas  scheme  on  a 
rectangular  domain  can  be  found  in  Turkel  and  Zwas  [1] 
and  Song  and  Tang  [13].  For  spherical  coordinated,  the 
CFL  condition  is  given  in  Navon  and  deVilliers  [5]: 

.  a  cos  9  pA\ 

At  * - - - . 

V^H  sin pkAX 

The  importance  of  this  explicit  scheme  is  also  in  paralleliza¬ 
tion.  It  is  always  easier  to  use  an  explicit  scheme  on  parallel 
computers  and  certainly  these  modifications  and  the  analy¬ 
sis  will  encourage  the  development  of  a  parallel  Turkel- 
Zwas  scheme. 

In  Section  2  we  present  the  modal  expansion  for  linear¬ 
ized  shallow  water  equations  in  spherical  coordinates  with 
no  mean  flow.  We  start  (Subsection  2.1)  by  obtaining  the 
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linearization  of  the  shallow  water  equations  on  the  sphere. 
Then  (Subsection  2.2)  we  describe  the  modal  expansion 
for  the  linearized  system  as  obtained  by  Longuet-Higgins 
[12].  In  Section  3  we  present  the  modal  expansion  for  the 
linearized  shallow  water  equations  in  spherical  coordinates 
discretized  using  the  Turkel-Zwas  finite  difference 
scheme.  We  start  this  section  by  describing  the  Turkel- 
Zwas  scheme  (Subsection  3.1).  The  modal  expansion  is 
given  in  Subsection  3.2.  As  a  result  of  this  analysis,  we 
conclude  that  certain  staggering  is  necessary.  A  modified 
staggered  method,  different  from  those  suggested  by  Song 
and  Tang  [13],  is  given  in  Section  4.  Again  this  section 
is  subdivided  into  two  subsections.  The  first  subsection 
describes  the  staggered  discretization  and  the  second  gives 
the  modal  expansion.  In  Section  5  we  compare  the  modal 
expansion  for  the  discrete  and  continuous  cases.  In  Section 
6  we  present  the  numerical  results  of  this  study  showing 
the  benefits  gained  by  staggering  the  Turkel-Zwas  scheme 
for  the  shallow  water  equations,  and  in  Section  7  we  give 
the  concluding  remarks. 

2.  SHALLOW  WATER  EQUATIONS 

2.1.  Linearization 

The  shallow  water  equations  in  spherical  coordinates 
are  given  by 


The  two-dimensional  linearized  shallow  water  equations 
(LSW)  with  no  mean  flow  in  spherical  coordinates  are 


dll  +  g  dh 
dt  a  cos  6  dX 


=  0 


dv 

dt 


+/u  +  £— =  o 

J  a  dO 


(5) 

(6) 


dh  |  H 
dt  a  cos  9 


du  d  . 

Ht  +  M(vCOSl)} 


=  0, 


(7) 


where  H  is  the  mean  height  of  the  surface  and  h  is  the 
perturbation.  In  the  next  subsection  we  describe  the  modal 
expansion  given  by  Longuet-Higgins  [12]. 

2.2.  Modal  Expansion  to  the  LSW  ( Longuet-Higgins ) 

For  the  case  of  interest  here,  /is  not  constant.  We  follow 
the  discussion  in  Longuet-Higgins  [12]  on  Laplace’s  tidal 
equations,  which  are  (5)-(7).  We  should  remark  here  that 
the  variable  6  in  Longuet-Higgins  [12]  is  trl2  —  9.  We  seek 
periodic  solutions  to  (5)-(7)  which  are  proportional  to 
ei(m\-ct ),  where  m  js  a  nonnegative  integer  and  c  a  constant 
nonzero  frequency.  Substituting  such  an  exponential  in 
(5)-(7),  we  get  im  instead  of  d/d  A,  and  —ic  instead  of  d/dt. 

Longuet-Higgins  [12]  introduced  functions  analogous  to 
the  velocity  potential  (<F)  and  stream  function  ('T),  such 
that 


dll  1 


dt  a  cos  6 


dll  „  dll 

u  —  +  v  cos  9 — - 
d\  d9 


-|/+- tan  y)tiH - =  0 

'a  a  cos  9  dX 


dv_  _ 1_ 

dt  a  cos  9 


dv  dv 

u - b  v  cos  9  — 

3A  d9 


u 

-1 

a 


dh 


1 


d  3 

zpr(hu)  +  ~(hv  cos  9) 
6X  d9 


dt  a  cos  9 
Here,  /is  the  Coriolis  parameter  given  by 


(1) 


-  1  + 

COS  9  dX  d9 

_  1  dW 

d9  COS  9  dX  ' 


(8) 

(9) 


Then  V2<b.  V2'!'  are  the  divergence  and  vorticity,  where 
V2  denotes  the  horizontal  Laplacian  operator, 


(10) 


L+S^-O 

/  ci  86 

(2) 

V2  =  1 

1^1 

<  ad\ 

cos  0  — 

1  d2  ’ 

cos  6 

dO' 

V  dOj 

COS  9 dX2 

=  0. 


(3)  Longuet-Higgins  has  shown  that 


-  V2d>  +  20  sin  9V2W  +  20  cos  9u  +  £  V2/;  =  0  (11) 
dt  a  v 


/  =  20  sin  9, 


(4) 


-  V2^  -  20  sin  9V2®  -  20  cos  9v  =  0.  (12) 

dt 


where  O  is  the  angular  speed  of  the  rotation  of  the  earth, 
h  is  the  height  of  the  homogeneous  atmosphere,  u  and  v  are 
the  zonal  and  meridional  wind  components,  respectively, 
9  and  A  are  the  latitudinal  and  longitudinal  directions, 
respectively,  a  is  the  radius  of  the  earth,  and  g  is  the  gravita¬ 
tional  constant. 


Equation  (7)  can  also  be  written 
dll  H 

—  +  —  V2<F  =  0.  (13) 

dt  a  v  ' 

We  substitute  for  u  and  v  from  (8)  and  (9)  to  get 
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-V2  +  20^-1$ 

dt  d\ ) 


+  20  (  sin  0V2  +  cos  6  — W  =  -  ^  V2h  (14) 

'  d\  a  v  ’ 


-V2  +  2a—  )  V 
dt  dX 


In  the  next  section,  we  follow  a  similar  derivation  for  the 
semidiscretization  given  by  the  Turkel-Zwas  scheme.  We 
will  show  how  Eqs.  (18)-(20)  and  (23)  will  be  affected  by 
the  discretization  and  prove  that  staggering  is  necessary. 
It  should  be  noted  that  the  idea  can  be  applied  to  any 
numerical  (semidiscrete)  scheme. 

3.  UNSTAGGERED  TURKEL-ZWAS  SCHEME 


2fl  (sin  OV2  +  cos  =  0.  (15) 


We  seek  solutions  to  these  equations  (13)— (15)  propor¬ 
tional  to  e'(mA-c0,  where  m  is  a  nonnegative  integer  and  c 
denotes  the  radian  frequency.  For  convenience  we  use  the 
operator  D  as 


3.1.  Discretization 

Given  a  constant  a,  0  <  a  <  1,  the  Turkel-Zwas  scheme 
for  the  nonlinear  shallow  water  equations  in  spherical  coor¬ 
dinates  takes  the  form 


u'kj  =  u'kj  -  <T 


cos  0, 


Oi+lJ  -  Mfc-lj)  +  vlkJ(ulkJ+1  -  l4,;-l) 


D -(1 V 


(16) 


+ 


g 


p  cos  6t 


( h'k+pj  ~  K-p.j) 


where 

f jl  =  sin  9. 

Then  Eqs.  (13)— (15)  become 


(AV2  —  m)  <&  +  (pV2  +  D)  iW  =  -  ^  V2/t  (18) 


ich  =  —  V2<F, 
a 


where  now 


v2  =  — 

dp. 


dp  J  1  —  p2 


and 


*  c 
~  2  a' 


Upon  eliminating  h  from  (18)  and  (20)  we  get 


+  2A  t 


u'kj  ( 


(l-«)  [  +  tan  6j  ui. 


(17) 


+  2  yi  +  ~^tan6nVk+pj 


+  2\fi  +  —<rtanei>Vk-p’1' 


(25) 


(AV2  -  m)  W  +  (pV2  +  D)  =  0  (19) 

(20) 


vk,j  —  vkJ  cr 


;  ( v'k+u  ~  v'k- 1,/)  +  As  (v'kj+ 1  -  v'kj. i) 


n  ^ lk-i+‘ 1  hkJ-q) 

(/ 


2A  t 


(21) 


ut 

(1  -  a)\^fj  +  -^tan  Ojj  u'kJ 
+  f  +  ^ptan  ulkJ+q 


(22) 


+  2  1 fhq  +  tan  °hq )  U‘k-i-q 


(26) 


hkj  =  hlkj  -  oT — 'hL-{ h'k+1j  -  hlk-i,j)  +  vkj  (hkJ+1  -  hlkj-i) 

l^COS  U; 


AV2  -  m  +  -A-  V4  <D  +  (pV2  +  D)  i<fr  =  0,  (23) 

eA 


+ 


hi 


fcj- 


COS  9. 


(1  -  a)(uk+pj  -  u'k-p.j) 


where 


T  2  ( tt k+p.j+q  tlk-p  j+q  T  lik+p  j-q  llk-p  j-c^) 


4  a2a2 
gH 


(24) 


+ 


hi 


Ly- 


COS  0; 


(1  -  a)(ui,/+9  cos  0j+q  -  vkJ-q  cos  0,_9) 
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+  2  {v'k+P,j+q  cos  6j+q  -  vi+pj-q  cos  0j-q 


+  vk-pJ+q  cos  0j+q  -  vk-pJ-q  cos  6hq) 


(27) 


Note  that  there  is  a  typographical  error  in  Eq.  (11a)  of 
Turkel  and  Zwas  [1],  which  is  our  Eq.  (29).  We  have  also 
modified  (to  get  a  symmetric  approximation,  as  suggested 
by  Neta  [9]  for  a  rectangular  domain)  the  right  hand  side 
of  (11c)  in  Turkel  and  Zwas  [1],  which  is  (31)  here. 


where 


A  t 


At 


a  AA  a  AO' 


(28) 


We  note  that  the  A,  6  components  of  the  pressure  gradi¬ 
ent  and  divergence  terms  are  differenced  over  points p  AX, 
q  Ad  (respectively)  away  rather  than  nearest  neighbors.  In 
this  way,  the  terms  associated  with  the  gravity  waves  (which 
carry  little  energy)  can  be  approximated  less  accurately 
and  thus  a  larger  time  step  can  be  allowed.  For  a  =  0, 
Turkel  and  Zwas  showed  that  there  is  a  deterioration  of 
accuracy  for  p  larger  than  one.  For  a  =  §  the  geostrophic 
balance  and  the  incompressibility  condition  are  satisfied 
to  a  higher  order  in  the  Cartesian  coordinate  case  (see 
Turkel  and  Zwas  [1],  Navon  and  de  Villiers  [5]  and  Navon 
and  Yu  [6]).  For  the  spherical  coordinates  case,  again  Tur¬ 
kel  and  Zwas  show  that  the  choice  a  =  5  was  necessary  to 
preserve  accuracy  when  p  or  q  is  greater  than  one. 

The  Turkel-Zwas  scheme  for  the  approximation  of  (5)- 
(7)  can  be  written  as 


du 

dt 


g 


k,j  2aAX  p  cos  0j 
+  fi 


(hk+p,j  hk-p,j) 


(1  -  a)vkJ  +  -  ( vk+pJ  +  vk-pj ) 


dv 

dt 


k,j  2a  A  0  q 


1  g 

(hk,j+q  h k . /  q ) 


3.2.  Modal  Expansion  for  the  Discretized  System 

We  follow  the  method  of  Longuet-Higgins  to  formulate 
the  discrete  eigenvalue  problem.  Now,  define  two  functions 
analogous  to  the  velocity  potential  and  stream  function 
such  that  at  each  point 


Ukj 


1 


cos  0, 


8P®k,  +  8,fVkj 


V  kj  8qQkj 


cos  Oj 


(32) 

(33) 


which  are  the  discrete  analog  of  (8)  and  (9),  where 

®k+p/2  j  ~  Q’k-pl 2  j 


Sfikj  = 


8q<S>kj  = 


p  AA 

(tV  j+q/2  ~  (lV  j-qi2 

q~A0  ' 


(34) 

(35) 


It  can  be  shown  that  V2d>  and  V2T  satisfy  the  discrete 
analog  of  (10);  i.e.. 


V2<IV,  =  — 82,<£>ki  +  — l-7r  Sq( cos  0y593>fcy),  (36) 


(29) 


(X 

(1  —  OL )  fjUjcj  +  —  ( fj+qUk,j+q  +  fj-qUk,j-q) 

(30) 

1 

z 

cos  Oj 

Oh 

dt 


1  H 


k,j  2 a  cos  Oj  L 


(l-«) 


p  AA 


k ’  cos2  Oj  p  ki  cos  Oj 


and  similarly  for  T. 

Now  combine  the  right  hand  sides  of  Eqs.  (29)  and  (30) 
in  the  following  way: 


(29)  A:  +  p/2  f  -  (29)k  -  p/2  j 
p  AA 

(30 ) kj  +  q/2  cos  0j+ql2  -  (30)kj.ql2  cos  0j.q/2 

qAO 


+ 


+ 


Vkj+q  cos  0j+q  -  v k  j-q  cos  0j-, 
qAO 


^  (X  (llk+pq+q  tlk-pj+q  +  ltk+pj-q 


Uk- 


+ 


+ 


2  \  p  AA  p  AA 

V k+p,j+q  COS  0j+q  V  k+p,j-q  COS  0j-q 

q~A0 

V k-pj+q  COS  Oj+q  V k-pj-q  COS  0j-q 

q ~A0 


(31) 


After  a  lengthy  algebraic  manipulation,  one  has  to  0((q 
AO)2),  which  agrees  with  (11), 

~  V2<t >kj  +  20  sin  OjV^kj  +  2D  cos  6 )ukj  +  -  S2hkj  =  0.  (37) 

Taking 


4  This  means  that  we  have  to  evaluate  (29)  at  the  point  k  +  p/2  j 
instead  of  k  j  and  similarly  for  the  other  terms  in  this  equation  and  the  next. 
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1 


cos  6j 


+ 


(30)k+p/2j-(30)k-p/2j 
p  AA 

(29) k  j  +  ql 2  cos  0j+q/2  -  (29)kj-q,2  cos  fl;-g/2 


qA9 

one  has  (note  the  similarity  to  (12)) 
d 


4.  STAGGERED  TURKEL-ZWAS  SCHEME 

4.1.  Discretization 

As  was  shown  in  the  previous  section,  it  is  necessary  to 
stagger  the  grid.  The  staggered  version  of  the  Turkel-Zwas 
scheme  takes  the  form 


dt 


V2xVkj  —  20  sin  6>;-V23>*;-  —  20  cos  6jvkj  =  0.  (38) 


ulkj  =  u'kj  ~(r 


u{. 


cos  9, 


—  ( Mfc+ij  -  wjt-ij)  +  v'kj(u'kj+1  -  m*j-i) 


Substituting  m*j,  vk,  from  (32)  and  (33)  into  (37)  and  (38) 
one  obtains  (similarly  to  (14)  and  (15)) 


2g 


—  V2  +  20, Sp  j  (h/c;  +  20  (sin  0,Y2 


+  cos  07<59)  +  £  V2/zfc7  =  0 


(39) 


+  p  cos  9,  (h'k-pl2’1  h'k~pl2’’\ 

+  2A t  (1  -  a)(^fj+1-^- tan  9^jvlkJ 

+  f  (  fj ■  +  tan  fy) 


V2  +  20<5p^  'P*;  —  20  (sin  0,V2  +  cos  9jSq)  Qkj  =  0.  (40) 

Equation  (31)  can  also  be  written  (to  second  order  in  q  A  9) 


+  «  tan^K-^ 


(43) 


j./+l  7  >  /  1  

vk,j  —  vkJ  cr 


cos  9, 


(Vk+lj  ~  Vlk-hj)  +  VlkJ(vlkJ+ 1  -  17^-i) 


dhkj  _ 

h  r  1 

dt 

a  (cos2  8j 

( [h/cJ+q/2  h/c,j-q/2) 


Slp8p<\>k  j 


+  co^ll  ^(cOS  W**') 


+  -  82q8p)^kj 


(41) 


2  At 


(1  -  a)(^fj  +  ^tan  0^  u‘kJ 


Note  that  (41)  is  slightly  different  from  (13).  The  first  two 
terms  would  combine  to  V23>*7  if  we  had  Sp,  8q  instead  of 
d2 p,  8lq,  respectively.  In  this  case,  the  third  term  would 
then  vanish.  It  may  be  easier  to  rewrite  (41)  as 


,  OL  /  UkJ+ql 2  \  1 

+  2  I  Jj+ql 2  +  ^  tan  0j+ql2  )  Uk,j+q/2 


.Oil  .  Ukj-q/2  I  ; 

+  2  1  fj-qn  +  - - - tan  0/-«/2  )  Uk,j-ql2 


(44) 


y  J  (cos  9j 


(h'k+i,j  ~  h'k- ij7) 


a/tq-  _  H  ,  1 


^  cqs2  q  ( ^2P  8p)  Sp&k 


+  Vlkj(h‘kj+ 1  -  /j*,;-l) 

(1  Oi)(llk+p!2,j  tlk-p/2J ) 


2hij_ 

COS  0, 


+  ^(S2‘i-S^cose’S^i)  (42) 


(  U/(.  1  p/2J  \  q!2  Uk  p/2,j+q  12 


2 - a  (82 p8q  ~  ^A)'^ j  f, 

cos  Oj  1 

where  the  first  term  is  now  the  same  as  (13)  and  the  others 
represent  the  perturbation.  These  other  terms  actually  sug¬ 
gest  that  one  should  take  Uk±P/2j,  vkj±qn  in  (31),  i.e.,  stag¬ 
gering  the  variable  h. 


nk+p/2,j-q/2  uk-pl2,j-qn) 


1 

J  p 


+  ■ 


2  hi 


h,i 


cos  0, 


(1  -  a)(v'kJ+q/ 2  cos  ej+qn 


"j  l 

—  V'k,j-qn  cos  9j-q/ 2) 
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( Vlk+p/2,j+ql2  COS  0j+q/2 


Vk+pl2,j-q/2  COS  Oj-q/2 


•  U 


+  V k-p/2,j+ql2  COS  dj+q/2  Vk-p/2,j-q/2  COS  0j-q/ 2) 

(45) 

where  cris  given  by  (28).  If  p  and  q  are  even  integers,  the 
staggered  grid  does  not  require  any  intermediate  values. 
In  Fig.  1,  we  show  the  location  of  the  grid  points  used 
by  the  unstaggered  Turkel-Zwas  scheme  in  each  of  the 
equations  (for  p  =  q  =  3).  For  the  staggered  grid  (again 
p  =  q  =  3)  intermediate  values  are  required  for  h  in  the 
momentum  equation  and  for  u  and  v  in  the  continuity 
equation.  It  is  interesting  to  see  the  difference  between 
the  staggered  grid  with  p  =  q  =  2  and  the  unstaggered 
grid  with  p  =  q  =  1.  One  may  feel  that  if  we  have  pi 2, 
ql 2  instead  of  p ,  q  respectively,  the  two  above  mentioned 
grids  are  the  same.  The  truth  is  that  in  the  momentum 
equation  only  h  moves  closer  to  the  center,  and  thus  the 
grids  are  different,  see  Fig.  2.  In  the  h  equation,  both  u 
and  v  move  closer  to  the  center  and  the  staggered  grid 
with  p  =  q  =  2  is  now  the  same  as  the  unstaggered  grid 
with  p  =  q  =  1,  see  Fig.  3.  Note  that  the  grid  for  the  h 
equation  is  identical  to  that  in  Fig.  2  and  thus  not  repeated. 


•  u,v,h»  •  u  •  u,v  •  u  • 


•  u 

u  equation 


•  u,h 


•  v 


•  V  •  U,V  •  V 


•  V 


4.2.  Modal  Expansion  for  the  Staggered  System 

We  now  rewrite  (29)-(31)  for  the  staggered  grid  (semi¬ 
discrete  staggered  Turkel-Zwas  scheme) 


•  u,h 

v  equation 


dll 

dt 


y  a  AS  p  cos  0j 


(hk+p/2,j  hk-p/2,j) 


+  fj 


a 


(1  -  a)vkJ  +  -  (vk+pJ  +  vk-pJ ) 


(46) 


dv 

dt 


k,j  a  A6q 


\  g 

—  (h-kj+q/2  —  hk,j-ql2 ) 


(1  a)fjl*k,j  T  2  (fj+qllk,j+q  T  fj-qUk,j-q) 


(47) 


dh 

dt  kj 


1  H 


a  cos  6; 


(1"«) 


uk+pl2,j  uk-pl2,j 

p  AA 


^  C k..j+q!2  COS  0j+q/2  t) kj- 

'  n  SO 


-„n  COS  0i . 


Q!  /  Uk+pl2,j+ql2  Uk-pl2,j+q/2 

2  \  p  AA 

Uk+p/2,j-q/2  ~  Uk-pl2,j-ql2 

p  AA 


•  u,v  •  •  •  v  •  • 


•  h 


•  u  •  •  h  •  u,v  •  h  • 


•  h 


•  U,V  •  •  •  V  •  • 

h  equation 

FIG.  1.  Unstaggered  grid  with  p  =  q  =  3. 


•  U,V,h 


•  U,V 


•  U 


•  U,V 
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+ 


+ 


V k+p/2,j+q/2  COS  Oj+q/ 2  V k+p/2,j-q/2  COS  Oj-q/ 2 

-qAO 

V k-p/2,j+q/2  COS  Oj+q/ 2  V k-p/2,j-q/2  COS  0j-q/2 

ql^O 


(AV  +  iSp )  (\>kj  +  (/uV  +  Vl  —  +  8q )  i'Vkj 


=  -i-^—V2In 
2  aCL  ki 


(50) 


(48) 


(AV2  +  iSp)  'Vkl  -  (+T  +  Vl^?Sq)  i$kj  =  0  (51) 


Equations  (39),  (40)  are  not  affected  (to  second  order 
in  q  AO)  by  the  staggering.  Equation  (42)  simplifies  to 


hki=  -  —  V2^,. 

7  ac  1 


(52) 


Eliminating  hkj  from  (50)  by  using  (52)  we  have 


dhkj  _  H, 
dt  a 


(49) 


(AV2  +  iSp  +  ^  V4)  <S>kj 


We  seek  solutions  to  (39),  (40),  and  (49)  proportional 
to  e~,c\  where  c  denotes  the  radian  frequency.  These  equa¬ 
tions  become 


+  (/rV2  +  Vl  -  ix28q)  Wkj  =  0  (53) 


(/rV2  +  i<bk]  -  (AV2  +  iSp)  Vkj  =  0.  (54) 

Suppose  <$kj,  xVkj  depend  on  \k  in  the  following  way. 


•  u 


•  u,v  •  u,h  •  u.v  •  u,h  •  u,v 


u  equation 


•  u 


•  u 


$a7-  =  FjeimXk 
+kj  =  G.eim+ 

where  m  is  a  nonnegative  integer;  then 

ft  \T/\  .  =  ft  Ci  <>lmXk  =  ^;+‘?/2 _ ^i-ql2  „unki. 

9  k>  q  >  qAO 

SpQkj  =  imf]Fje,mXk 


(55) 

(56) 


.  girriAfc 


(57) 

(58) 


•  v. 


and 


•  v  •  u,v  •  v 


•  v,h 


•  u 


v  equation 

•  u.v  •  v,h  •  u.v 

•  u,h  •  u,v  •  u,h 


•  u.v  •  v,h  •  u,v 
h  equation 


•  u 

•  u.v,h*  u,v  •  u,v,h 

•  u 

u  equation 

•  u,v,h 

•  v  •  u,v  •  v 

•  u,v,h 

v  equation 


FIG.  2.  Staggered  grid  with  p  =  q  =  2. 


FIG.  3.  Unstaggered  grid  with  p  =  q  =  1. 
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V<u7  = 


-(my])2  |  1 

_  1  -  V  Vl  -  p? 


Sq(V]  -  ]x2Sq) 


FfiimXk. 


(59) 


where 


.  _  sin(mp  AA/2) 

^  ( mp  AA/2) 

Upon  substituting  (58)  in  (53)  and  (54),  we  get  the  semidis¬ 
crete  analogs  of  (23)  and  (19),  respectively. 

5.  COMPARISON  OF  DISCRETE  AND  CONTINUOUS 
MODAL  EXPANSIONS 

In  this  section,  we  compare  the  equations  used  for  the 
modal  expansion  in  the  continuous  case,  i.e.,  (23)  and  (19), 
with  the  corresponding  ones  for  the  Turkel-Zwas  discreti¬ 
zation.  These  equations  for  the  staggered  Turkel-Zwas 
scheme  are  (53)  and  (54),  respectively.  In  the  unstaggered 
case,  we  were  unable  to  obtain  such  a  system,  since  (41) 
is  different  in  form  from  (13).  This  difference  led  us  to  the 
development  of  the  staggered  scheme. 

Note  the  similarity  between  (53)  and  (23).  The  opera¬ 
tor  V2,  as  defined  by  (21),  in  (23)  is  replaced  by  its  discrete 
analog  defined  in  (36).  The  operator  Z),  defined  by  (16), 
is  now  discretized  by  Vl  —  p}Sq .  The  discrete  operator  iSp 
becomes  —  m  in  the  continuous  case.  The  same  analogy 
holds  between  (19)  and  (54). 


nr 


4>  =  5.768  X  104-^V 


sec~ 


u„  =  20 


m 


sec 

a  =  6.370  X  106  m 

H  =  7.292  X  HU5—. 

sec 


The  / 2  error  norms  for  the  height  and  velocity  are  computed 
in  the  manner 


,,  2  =  V/[(/r(A,  fl)-ftf(A,  ft))2] 
11 v/[(/a  a,  0))2] 


and 


V/[(u(A,  9)  -  ue( A,  8))2  +  (u(A,  9)  -  ue(A,  ft))2] 
V/[(Me( A,  ft))2  +  (ue(A,  0))2] 


where 


1  C2.TI  jr/2 

^[/(A,  e)]  =  —  /(A,  0)  cos  0ddd\. 

477  J  0  J  -jf/2 


In  these  relations  the  variables  he ,  ue,  and  ve  are  considered 
to  be  the  exact  solution.  This  solution  is  obtained  by  a 
centered  in  space  and  time  algorithm  (a  leapfrog  scheme) 
on  a  128  X  64  grid  that  is  time  integrated  using  a  time  step 


6.  NUMERICAL  RESULTS 


In  this  section  we  present  the  numerical  results  for  the 
unstaggered  and  staggered  versions  of  the  Turkel-Zwas 
scheme.  The  initial  conditions  are  the  same  as  those  used 
by  McDonald  and  Bates  [14],  where  the  height  is  defined  as 


/? (A,  6, 0) 


\  _ 

-(<£>  +  2f lau0  sin3  6  cos  9  sin  A) 

o 


and  the  velocities  are  obtained  geostrophically  to  yield 

u( A,  9 , 0)  =  —3 ua  sin  9  cos2  0sin  A  +  u0  sin3  9  sin  A 
v  ( A,  9, 0)  =  u0  sin2  9  cos  A, 


where 


g  =  9.8 


m 

sec2 


FIG.  4.  The  height  at  the  initial  state  (t  =  0  h).  The  grid  is  128  X  64. 
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FIG.  5.  The  height  at  the  final  state  (t  =  24  h).  The  grid  is  128  X  64 
and  the  time  step  is  A t  =  15  s. 


of  15  s.  The  height  at  the  initial  (t  =  0  h)  and  final  ( t  = 
24  h)  states  is  shown  in  Figs.  4  and  5. 

The  percentage  change  in  total  available  energy  is  de¬ 
fined  as 


A  E 


I [E final]  ~  I[E  initial} 

initial] 


X  100%, 


where 


E( A,  9) 


h{ A,  9){u{ A,  0)2 


+  v{\,  9)2)  + 


U( A,  9) 


The  shallow  water  equations  are  integrated  in  time  for 
a  period  of  24  h  on  a  64  X  32  grid  on  a  Sparc  10.  Table  1 
shows  the  results  for  various  configurations  of  p  and  q 
using  different  time  steps  for  the  unstaggered  and  stag¬ 
gered  versions  of  the  Turkel-Zwas  scheme.  Recall  that  p 
and  q  refer  to  the  extension  of  the  differencing  stencil  in  the 
longitudinal  (A)  and  latitudinal  ( 9)  directions,  respectively, 
and  a  is  the  Pade-type  differencing  weight.  We  have  experi¬ 
mented  with  various  values  of  a,  namely  a  =  0,  J,  J,  §, 
I,  1.  We  have  found  that  a  =  J  is  the  best  choice.  The  first 
row  in  Table  1  represents  the  leapfrog  scheme,  which  only 
allows  a  time  step  of  100  s  due  to  the  very  restrictive 
CFL  condition  which  limits  all  explicit  algorithms.  For  the 
remainder  of  the  cases,  p  is  always  greater  than  q  because 
the  time  restriction  is  dictated  by  the  distance  between 
mesh  points  in  the  A  direction.  As  p  increases  so  too  does 
the  maximum  allowable  time  step  because  the  CFL  restric¬ 
tion  is  relaxed  due  to  the  coarser  differencing  stencil  used 
for  the  gravity  wave  terms.  For  typical  meteorological  con¬ 
ditions  the  inertial  gravity  wave  speeds  are  larger  than  the 
wind  velocities  while  most  of  the  energy  is  carried  by  the 
wind  velocities.  Thus  it  makes  sense  to  treat  the  gravity 
terms  less  accurately  by  using  a  coarser  grid  for  these  terms 
(see  Turkel  and  Zwas  [1]). 


TABLE  1 

The  l 2  Error  Norms  for  the  Height  h  and  the  Velocity  Field  u  for  the  Unstaggered  and  Staggered  Versions 
of  the  Turkel-Zwas  Scheme  for  a  64  X  32  Grid  after  24  h 


Staggered  Grid 

P 

<7 

a 

At 

(s) 

WV 

(X  10~4) 

ll“llr 

(x  ur3) 

A  E 

(%) 

CPU 

(s) 

No 

1 

1 

0 

100 

1.177 

3.722 

-0.09 

240 

No 

2 

1 

1/3 

200 

1.187 

3.830 

-0.06 

119 

No 

3 

2 

1/3 

200 

2.367 

8.018 

-0.15 

119 

Yes 

3 

2 

1/3 

200 

1.221 

4.096 

-0.03 

119 

No 

4 

2 

1/3 

200 

2.406 

8.387 

-0.15 

119 

Yes 

4 

2 

1/3 

200 

1.269 

4.478 

0 

119 

No 

3 

1 

1/3 

300 

1.193 

3.931 

-0.06 

83 

No 

5 

2 

1/3 

300 

2.475 

8.844 

-0.15 

83 

Yes 

5 

2 

1/3 

300 

1.344 

5.081 

+0.03 

83 

No 

6 

2 

1/3 

300 

2.592 

9.471 

-0.15 

83 

Yes 

6 

2 

1/3 

300 

1.467 

5.977 

+0.06 

83 

No 

4 

1 

1/3 

400 

1.211 

4.149 

-0.06 

61 

No 

7 

2 

1/3 

400 

2.735 

10.06 

-0.15 

61 

Yes 

7 

2 

1/3 

400 

1.634 

7.131 

+0.09 

61 

No 

8 

2 

1/3 

400 

2.881 

10.90 

-0.15 

61 

Yes 

8 

2 

1/3 

400 

1.871 

8.585 

+0.15 

61 

Note.  The  parameter  A  E  represents  the  percentage  change  of  available  energy  of  the  system. 
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Table  1  shows  that  for  lower  values  of  p  and  q ,  say 
(p  =  2,  q  =  1)  for  the  unstaggered  and  (p  =  3,  q  =  2)  for 
the  staggered  and  At  =  200,  the  unstaggered  and  staggered 
cases  yield  comparable  error  norms.  As  p  and  q  are  in¬ 
creased  beyond  these  values,  the  errors  increase  dramati¬ 
cally  for  the  unstaggered  case;  the  errors  for  the  staggered 
case  are  only  half  those  of  the  unstaggered  case.  Note  that 
for  the  same  time  step,  the  staggered  case  takes  no  more 
CPU  time  than  the  unstaggered  case.  In  addition,  the  stag¬ 
gered  case  where  (p  =  4,  q  =  2)  conserves  the  total  avail¬ 
able  energy.  As  p  increases,  the  staggered  case  no  longer 
conserves  the  available  energy  but  it  is  still  more  accurate 
than  the  unstaggered  case  for  the  same  values  of  p  and  q. 
Therefore  for  the  staggered  case  the  optimal  values  of  p 
and  q  must  lie  somewhere  in  the  vicinity  of  (p  =  3,  q  = 
2)  and  (p  =  4,  q  =  2).  By  comparing  the  error  norms  for 
the  unstaggered  case  for  the  values  (p  =  3,  q  =  1)  and 
(p  =  3,  q  =  2),  and  (p  =  4,  q  =  1)  and  (p  =  4,  q  =  2), 
we  see  that  the  increase  in  q  from  1  to  2  has  adverse  effects 
on  the  solution  accuracy.  Therefore  it  stands  to  reason  that 
for  the  staggered  case  a  better  solution  can  be  achieved 
by  the  values  (p  =  3,  q  =  1)  or  (p  =  4,  q  =  1).  The  difficulty 
with  this  case  is  that  now  interpolation  is  required  because 
the  staggering  for  odd  values  of  p  and  q  results  in  interme¬ 
diate  values.  Intermediate  values  in  the  longitudinal  direc¬ 
tion  pose  no  difficulty  and  in  fact  are  handled  by  linear 
interpolation  in  this  paper.  On  the  other  hand,  intermedi¬ 
ate  values  in  the  latitudinal  direction  cause  problems  be¬ 
cause  interpolation  is  now  required  at  the  poles  where  the 
wind  velocities  are  undefined.  This  case  is  not  included  in 
this  paper. 

For  completeness  let  us  review  the  case  Af  =  200  once 
again.  The  staggered  case  (p  =  4,  q  =  2)  yields  better 
results  than  the  unstaggered  case  (p  =  4,  q  =  2).  However, 
the  unstaggered  case  allows  a  time  step  of  400  s.  The  stag¬ 
gered  case,  on  the  other  hand,  cannot  use  so  large  a  time 
step  because  the  CFL  condition  is  dictated  by  terms  using 
the  differencing  stencil  points  (p/2,  q/2)  instead  of  (p,  q). 
This  means  that  there  is  a  trade-off  between  computing 
time  and  accuracy;  for  a  given  value  of  p  and  q,  the  unstag¬ 
gered  case  allows  a  larger  time  step  than  the  staggered 
case  thereby  requiring  less  CPU  time  at  the  cost  of  lower 
accuracy.  To  clarify  the  point,  we  recall  that  Song  and 
Tang  have  proved  that  basically  the  stability  condition  of 
the  unstaggered  case  is 

Af  _  p 

d  V2 ~gh 

Since  in  our  case,  we  have  the  points  half  way  the  original 
(unstaggered)  distance,  we  replace  the  numerator  by  p/2, 
which  is  the  same  as  halving  the  value  of  Atld.  For  the 
spherical  case,  a  similar  dependence  was  proved  by  Turkel 
and  Zwas  [1]  (see  our  Introduction). 


7.  SUMMARY  AND  CONCLUSIONS 

The  analysis  of  the  spherical  coordinates  shallow  water 
equations  model  shows  that  the  T-Z  scheme  must  be  stag¬ 
gered  to  get  eigenvalues  and  eigenfunctions  approaching 
those  of  the  continuous  case.  The  importance  of  such  an 
analysis  is  the  fact  that  it  is  valid  for  nonconstant  coeffi¬ 
cients  and  thereby  applicable  to  any  numerical  scheme.  In 
addition,  numerical  experiments  are  conducted  illustrating 
the  benefits  of  staggering  the  original  scheme.  The  numeri¬ 
cal  experiments  show  that  for  given  values  of  p  and  q  the 
staggered  case  performs  better  than  the  unstaggered  case. 
However,  the  staggered  case  requires  half  the  time  step 
of  the  unstaggered  case.  This  means  that  there  is  a  trade¬ 
off  between  efficiency  and  accuracy;  for  a  given  value  of 
p  and  q,  the  unstaggered  case  allows  a  larger  time  step 
than  the  staggered  case  thereby  requiring  less  CPU  time 
at  the  cost  of  lower  accuracy.  These  experiments  also  show 
that  the  best  results  for  the  staggered  case  are  obtained 
with  the  values  (p  =  3,  q  =  2)  and  {p  =  4,  q  =  2).  Further¬ 
more,  the  experiments  suggest  that  better  results  may  be 
obtained  by  the  configurations  (p  =  3,  q  =  1)  and  (p  = 
4,  q  =  1)  but  interpolation  is  required  in  the  latitudinal 
direction  because  odd  values  of  p  or  q  result  in  differencing 
points  that  do  not  lie  on  grid  points  (intermediate  values). 
Interpolation  in  the  longitudinal  direction  is  straightfor¬ 
ward  and  is  handled  by  linear  interpolation  in  this  paper, 
while  in  the  latitudinal  direction  it  is  no  longer  trivial  be¬ 
cause  intermediate  values  may  fall  on  the  poles  where  the 
wind  velocities  are  undefined. 

Software  for  both  the  unstaggered  and  staggered  Tur- 
kel-Zwas  schemes  for  the  approximation  of  the  shallow 
water  equations  in  spherical  coordinates  is  available  at 
URL  http://math.nps.navy.mil/~bneta. 
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